1//////////////////////////////////////////////////////////////////////
   2// LibFile: paths.scad
   3//   A `path` is a list of points of the same dimensions, usually 2D or 3D, that can
   4//   be connected together to form a sequence of line segments or a polygon.
   5//   A `region` is a list of paths that represent polygons, and the functions
   6//   in this file work on paths and also 1-regions, which are regions
   7//   that include exactly one path.  When you pass a 1-region to a function, the default
   8//   value for `closed` is always `true` because regions represent polygons.  
   9//   Capabilities include computing length of paths, computing
  10//   path tangents and normals, resampling of paths, and cutting paths up into smaller paths.  
  11// Includes:
  12//   include <BOSL2/std.scad>
  13// FileGroup: Advanced Modeling
  14// FileSummary: Operations on paths: length, resampling, tangents, splitting into subpaths
  15// FileFootnotes: STD=Included in std.scad
  16//////////////////////////////////////////////////////////////////////
  17
  18// Section: Utility Functions
  19
  20// Function: is_path()
  21// Synopsis: Returns True if 'list' is a path.
  22// Topics: Paths
  23// See Also: is_region(), is_vnf()
  24// Usage:
  25//   is_path(list, [dim], [fast])
  26// Description:
  27//   Returns true if `list` is a path.  A path is a list of two or more numeric vectors (AKA points).
  28//   All vectors must of the same size, and may only contain numbers that are not inf or nan.
  29//   By default the vectors in a path must be 2d or 3d.  Set the `dim` parameter to specify a list
  30//   of allowed dimensions, or set it to `undef` to allow any dimension.  (Note that this function
  31//   returns `false` on 1-regions.)  
  32// Example:
  33//   bool1 = is_path([[3,4],[5,6]]);    // Returns true
  34//   bool2 = is_path([[3,4]]);          // Returns false
  35//   bool3 = is_path([[3,4],[4,5]],2);  // Returns true
  36//   bool4 = is_path([[3,4,3],[5,4,5]],2);  // Returns false
  37//   bool5 = is_path([[3,4,3],[5,4,5]],2);  // Returns false
  38//   bool6 = is_path([[3,4,5],undef,[4,5,6]]);  // Returns false
  39//   bool7 = is_path([[3,5],[undef,undef],[4,5]]);  // Returns false
  40//   bool8 = is_path([[3,4],[5,6],[5,3]]);     // Returns true
  41//   bool9 = is_path([3,4,5,6,7,8]);           // Returns false
  42//   bool10 = is_path([[3,4],[5,6]], dim=[2,3]);// Returns true
  43//   bool11 = is_path([[3,4],[5,6]], dim=[1,3]);// Returns false
  44//   bool12 = is_path([[3,4],"hello"], fast=true); // Returns true
  45//   bool13 = is_path([[3,4],[3,4,5]]);            // Returns false
  46//   bool14 = is_path([[1,2,3,4],[2,3,4,5]]);      // Returns false
  47//   bool15 = is_path([[1,2,3,4],[2,3,4,5]],undef);// Returns true
  48// Arguments:
  49//   list = list to check
  50//   dim = list of allowed dimensions of the vectors in the path.  Default: [2,3]
  51//   fast = set to true for fast check that only looks at first entry.  Default: false
  52function is_path(list, dim=[2,3], fast=false) =
  53    fast
  54    ?   is_list(list) && is_vector(list[0]) 
  55    :   is_matrix(list) 
  56        && len(list)>1 
  57        && len(list[0])>0
  58        && (is_undef(dim) || in_list(len(list[0]), force_list(dim)));
  59
  60// Function: is_1region()
  61// Synopsis: Returns true if path is a region with one component.
  62// Topics: Paths, Regions
  63// See Also: force_path()
  64// Usage:
  65//   bool = is_1region(path, [name])
  66// Description:
  67//   If `path` is a region with one component (a 1-region) then return true.  If path is a region with more components
  68//   then display an error message about the parameter `name` requiring a path or a single component region.  If the input
  69//   is not a region then return false.  This function helps path functions accept 1-regions.
  70// Arguments:
  71//   path = input to process
  72//   name = name of parameter to use in error message.  Default: "path"
  73function is_1region(path, name="path") = 
  74     !is_region(path)? false
  75    :assert(len(path)==1,str("Parameter \"",name,"\" must be a path or singleton region, but is a multicomponent region"))
  76     true;
  77
  78
  79// Function: force_path()
  80// Synopsis: Checks that path is a region with one component.
  81// SynTags: Path
  82// Topics: Paths, Regions
  83// See Also: is_1region()
  84// Usage:
  85//   outpath = force_path(path, [name])
  86// Description:
  87//   If `path` is a region with one component (a 1-region) then returns that component as a path.
  88//   If path is a region with more components then displays an error message about the parameter
  89//   `name` requiring a path or a single component region.  If the input is not a region then
  90//   returns the input without any checks.  This function helps path functions accept 1-regions.
  91// Arguments:
  92//   path = input to process
  93//   name = name of parameter to use in error message.  Default: "path"
  94function force_path(path, name="path") =
  95   is_region(path) ?
  96       assert(len(path)==1, str("Parameter \"",name,"\" must be a path or singleton region, but is a multicomponent region"))
  97       path[0]
  98   : path;
  99
 100
 101/// Internal Function: _path_select()
 102/// Usage:
 103///   _path_select(path,s1,u1,s2,u2,[closed]):
 104/// Description:
 105///   Returns a portion of a path, from between the `u1` part of segment `s1`, to the `u2` part of
 106///   segment `s2`.  Both `u1` and `u2` are values between 0.0 and 1.0, inclusive, where 0 is the start
 107///   of the segment, and 1 is the end.  Both `s1` and `s2` are integers, where 0 is the first segment.
 108/// Arguments:
 109///   path = The path to get a section of.
 110///   s1 = The number of the starting segment.
 111///   u1 = The proportion along the starting segment, between 0.0 and 1.0, inclusive.
 112///   s2 = The number of the ending segment.
 113///   u2 = The proportion along the ending segment, between 0.0 and 1.0, inclusive.
 114///   closed = If true, treat path as a closed polygon.
 115function _path_select(path, s1, u1, s2, u2, closed=false) =
 116    let(
 117        lp = len(path),
 118        l = lp-(closed?0:1),
 119        u1 = s1<0? 0 : s1>l? 1 : u1,
 120        u2 = s2<0? 0 : s2>l? 1 : u2,
 121        s1 = constrain(s1,0,l),
 122        s2 = constrain(s2,0,l),
 123        pathout = concat(
 124            (s1<l && u1<1)? [lerp(path[s1],path[(s1+1)%lp],u1)] : [],
 125            [for (i=[s1+1:1:s2]) path[i]],
 126            (s2<l && u2>0)? [lerp(path[s2],path[(s2+1)%lp],u2)] : []
 127        )
 128    ) pathout;
 129
 130
 131
 132// Function: path_merge_collinear()
 133// Synopsis: Removes unnecessary points from a path.
 134// SynTags: Path
 135// Topics: Paths, Regions
 136// Description:
 137//   Takes a path and removes unnecessary sequential collinear points.  Note that when `closed=true` either of the path
 138//   endpoints may be removed.  
 139// Usage:
 140//   path_merge_collinear(path, [eps])
 141// Arguments:
 142//   path = A path of any dimension or a 1-region
 143//   closed = treat as closed polygon.  Default: false
 144//   eps = Largest positional variance allowed.  Default: `EPSILON` (1-e9)
 145function path_merge_collinear(path, closed, eps=EPSILON) =
 146    is_1region(path) ? path_merge_collinear(path[0], default(closed,true), eps) :
 147    let(closed=default(closed,false))
 148    assert(is_bool(closed))
 149    assert( is_path(path), "Invalid path in path_merge_collinear." )
 150    assert( is_undef(eps) || (is_finite(eps) && (eps>=0) ), "Invalid tolerance." )
 151    len(path)<=2 ? path :
 152    [
 153      if(!closed) path[0],
 154      for(triple=triplet(path,wrap=closed))
 155        if (!is_collinear(triple,eps=eps)) triple[1],
 156      if(!closed) last(path)
 157    ];
 158
 159
 160// Section: Path length calculation
 161
 162
 163// Function: path_length()
 164// Synopsis: Returns the path length.
 165// Topics: Paths
 166// See Also: path_segment_lengths(), path_length_fractions()
 167// Usage:
 168//   path_length(path,[closed])
 169// Description:
 170//   Returns the length of the path.
 171// Arguments:
 172//   path = Path of any dimension or 1-region. 
 173//   closed = true if the path is closed.  Default: false
 174// Example:
 175//   path = [[0,0], [5,35], [60,-25], [80,0]];
 176//   echo(path_length(path));
 177function path_length(path,closed) =
 178    is_1region(path) ? path_length(path[0], default(closed,true)) :
 179    assert(is_path(path), "Invalid path in path_length")
 180    let(closed=default(closed,false))
 181    assert(is_bool(closed))
 182    len(path)<2? 0 :
 183    sum([for (i = [0:1:len(path)-2]) norm(path[i+1]-path[i])])+(closed?norm(path[len(path)-1]-path[0]):0);
 184
 185
 186// Function: path_segment_lengths()
 187// Synopsis: Returns a list of the lengths of segments in a path.
 188// Topics: Paths
 189// See Also: path_length(), path_length_fractions()
 190// Usage:
 191//   path_segment_lengths(path,[closed])
 192// Description:
 193//   Returns list of the length of each segment in a path
 194// Arguments:
 195//   path = path in any dimension or 1-region
 196//   closed = true if the path is closed.  Default: false
 197function path_segment_lengths(path, closed) =
 198    is_1region(path) ? path_segment_lengths(path[0], default(closed,true)) :
 199    let(closed=default(closed,false))
 200    assert(is_path(path),"Invalid path in path_segment_lengths.")
 201    assert(is_bool(closed))
 202    [
 203        for (i=[0:1:len(path)-2]) norm(path[i+1]-path[i]),
 204        if (closed) norm(path[0]-last(path))
 205    ]; 
 206
 207
 208// Function: path_length_fractions()
 209// Synopsis: Returns the fractional distance of each point along the length of a path.
 210// Topics: Paths
 211// See Also: path_length(), path_segment_lengths()
 212// Usage:
 213//   fracs = path_length_fractions(path, [closed]);
 214// Description:
 215//    Returns the distance fraction of each point in the path along the path, so the first
 216//    point is zero and the final point is 1.  If the path is closed the length of the output
 217//    will have one extra point because of the final connecting segment that connects the last
 218//    point of the path to the first point.
 219// Arguments:
 220//    path = path in any dimension or a 1-region
 221//    closed = set to true if path is closed.  Default: false
 222function path_length_fractions(path, closed) =
 223    is_1region(path) ? path_length_fractions(path[0], default(closed,true)):
 224    let(closed=default(closed, false))
 225    assert(is_path(path))
 226    assert(is_bool(closed))
 227    let(
 228        lengths = [
 229            0,
 230            each path_segment_lengths(path,closed)
 231        ],
 232        partial_len = cumsum(lengths),
 233        total_len = last(partial_len)
 234    )
 235    partial_len / total_len;
 236
 237
 238
 239/// Internal Function: _path_self_intersections()
 240/// Usage:
 241///   isects = _path_self_intersections(path, [closed], [eps]);
 242/// Description:
 243///   Locates all self intersection points of the given path.  Returns a list of intersections, where
 244///   each intersection is a list like [POINT, SEGNUM1, PROPORTION1, SEGNUM2, PROPORTION2] where
 245///   POINT is the coordinates of the intersection point, SEGNUMs are the integer indices of the
 246///   intersecting segments along the path, and the PROPORTIONS are the 0.0 to 1.0 proportions
 247///   of how far along those segments they intersect at.  A proportion of 0.0 indicates the start
 248///   of the segment, and a proportion of 1.0 indicates the end of the segment.
 249///   .
 250///   Note that this function does not return self-intersecting segments, only the points
 251///   where non-parallel segments intersect.  
 252/// Arguments:
 253///   path = The path to find self intersections of.
 254///   closed = If true, treat path like a closed polygon.  Default: true
 255///   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
 256/// Example(2D):
 257///   path = [
 258///       [-100,100], [0,-50], [100,100], [100,-100], [0,50], [-100,-100]
 259///   ];
 260///   isects = _path_self_intersections(path, closed=true);
 261///   // isects == [[[-33.3333, 0], 0, 0.666667, 4, 0.333333], [[33.3333, 0], 1, 0.333333, 3, 0.666667]]
 262///   stroke(path, closed=true, width=1);
 263///   for (isect=isects) translate(isect[0]) color("blue") sphere(d=10);
 264function _path_self_intersections(path, closed=true, eps=EPSILON) =
 265    let(
 266        path = closed ? list_wrap(path,eps=eps) : path,
 267        plen = len(path)
 268    )
 269    [ for (i = [0:1:plen-3]) let(
 270          a1 = path[i],
 271          a2 = path[i+1], 
 272          seg_normal = unit([-(a2-a1).y, (a2-a1).x],[0,0]),
 273          vals = path*seg_normal,
 274          ref  = a1*seg_normal,
 275            // The value of vals[j]-ref is positive if vertex j is one one side of the
 276            // line [a1,a2] and negative on the other side. Only a segment with opposite
 277            // signs at its two vertices can have an intersection with segment
 278            // [a1,a2]. The variable signals is zero when abs(vals[j]-ref) is less than
 279            // eps and the sign of vals[j]-ref otherwise.  
 280          signals = [for(j=[i+2:1:plen-(i==0 && closed? 2: 1)]) 
 281                        abs(vals[j]-ref) <  eps ? 0 : sign(vals[j]-ref) ]
 282        )
 283        if(max(signals)>=0 && min(signals)<=0 ) // some remaining edge intersects line [a1,a2]
 284        for(j=[i+2:1:plen-(i==0 && closed? 3: 2)])
 285            if( signals[j-i-2]*signals[j-i-1]<=0 ) let( // segm [b1,b2] intersects line [a1,a2]
 286                b1 = path[j],
 287                b2 = path[j+1],
 288                isect = _general_line_intersection([a1,a2],[b1,b2],eps=eps) 
 289            )
 290            if (isect 
 291                && isect[1]>=-eps
 292                && isect[1]<= 1+eps
 293                && isect[2]>= -eps 
 294                && isect[2]<= 1+eps)
 295                [isect[0], i, isect[1], j, isect[2]]
 296    ];
 297
 298// Section: Resampling - changing the number of points in a path
 299
 300
 301// Input `data` is a list that sums to an integer. 
 302// Returns rounded version of input data so that every 
 303// entry is rounded to an integer and the sum is the same as
 304// that of the input.  Works by rounding an entry in the list
 305// and passing the rounding error forward to the next entry.
 306// This will generally distribute the error in a uniform manner. 
 307function _sum_preserving_round(data, index=0) =
 308    index == len(data)-1 ? list_set(data, len(data)-1, round(data[len(data)-1])) :
 309    let(
 310        newval = round(data[index]),
 311        error = newval - data[index]
 312    ) _sum_preserving_round(
 313        list_set(data, [index,index+1], [newval, data[index+1]-error]),
 314        index+1
 315    );
 316
 317
 318// Function: subdivide_path()
 319// Synopsis: Subdivides a path to produce a more finely sampled path.
 320// SynTags: Path
 321// Topics: Paths, Path Subdivision
 322// See Also: subdivide_and_slice(), resample_path(), jittered_poly()
 323// Usage:
 324//   newpath = subdivide_path(path, n|refine=|maxlen=, [method=], [closed=], [exact=]);
 325// Description:
 326//   Takes a path as input (closed or open) and subdivides the path to produce a more
 327//   finely sampled path.  You control the subdivision process by using the `maxlen` arg
 328//   to specify a maximum segment length, or by specifying `n` or `refine`, which request
 329//   a certain point count in the output.
 330//   .
 331//   You can specify the point count using the `n` option, where
 332//   you give the number of points you want in the output, or you can use
 333//   the `refine` option, where you specify a resampling factor.  If `refine=3` then
 334//   the number of points would increase by a factor of three, so a four point square would
 335//   have 12 points after subdivision.  With point-count subdivision, the new points can be distributed
 336//   proportional to length (`method="length"`), which is the default, or they can be divided up evenly among all the path segments
 337//   (`method="segment"`).  If the extra points don't fit evenly on the path then the
 338//   algorithm attempts to distribute them as uniformly as possible, but the result may be uneven.
 339//   The `exact` option, which is true by default, requires that the final point count is
 340//   exactly as requested.  For example, if you subdivide a four point square and request `n=13` then one edge will have
 341//   an extra point compared to the others.  
 342//   If you set `exact=false` then the
 343//   algorithm will favor uniformity and the output path may have a different number of
 344//   points than you requested, but the sampling will be uniform.   In our example of the
 345//   square with `n=13`, you will get only 12 points output, with the same number of points on each edge.
 346//   .
 347//   The points are always distributed uniformly on each segment.  The `method="length"` option does
 348//   means that the number of points on a segment is based on its length, but the points are still
 349//   distributed uniformly on each segment, independent of the other segments.  
 350//   With the `"segment"` method you can also give `n` as a vector of counts.  This 
 351//   specifies the desired point count on each segment: with vector valued `n` the `subdivide_path`
 352//   function places `n[i]-1` points on segment `i`.  The reason for the -1 is to avoid
 353//   double counting the endpoints, which are shared by pairs of segments, so that for
 354//   a closed polygon the total number of points will be sum(n).  Note that with an open
 355//   path there is an extra point at the end, so the number of points will be sum(n)+1.
 356//   .
 357//   If you use the `maxlen` option then you specify the maximum length segment allowed in the output.
 358//   Each segment is subdivided into the largest number of segments meeting your requirement.  As above,
 359//   the sampling is uniform on each segment, independent of the other segments.  With the `maxlen` option
 360//   you cannot specify `method` or `exact`.    
 361// Arguments:
 362//   path = path in any dimension or a 1-region
 363//   n = scalar total number of points desired or with `method="segment"` can be a vector requesting `n[i]-1` new points added to segment i.
 364//   ---
 365//   refine = increase total number of points by this factor (Specify only one of n, refine and maxlen)
 366//   maxlen = maximum length segment in the output (Specify only one of n, refine and maxlen)
 367//   closed = set to false if the path is open.  Default: True
 368//   exact = if true return exactly the requested number of points, possibly sacrificing uniformity.  If false, return uniform point sample that may not match the number of points requested.  (Not allowed with maxlen.) Default: true
 369//   method = One of `"length"` or `"segment"`.  If `"length"`, adds vertices in proportion to segment length, so short segments get fewer points.  If `"segment"`, add points evenly among the segments, so all segments get the same number of points.  (Not allowed with maxlen.) Default: `"length"`
 370// Example(2D):
 371//   mypath = subdivide_path(square([2,2],center=true), 12);
 372//   move_copies(mypath)circle(r=.1,$fn=32);
 373// Example(2D):
 374//   mypath = subdivide_path(square([8,2],center=true), 12);
 375//   move_copies(mypath)circle(r=.2,$fn=32);
 376// Example(2D):
 377//   mypath = subdivide_path(square([8,2],center=true), 12, method="segment");
 378//   move_copies(mypath)circle(r=.2,$fn=32);
 379// Example(2D):
 380//   mypath = subdivide_path(square([2,2],center=true), 17, closed=false);
 381//   move_copies(mypath)circle(r=.1,$fn=32);
 382// Example(2D): Specifying different numbers of points on each segment
 383//   mypath = subdivide_path(hexagon(side=2), [2,3,4,5,6,7], method="segment");
 384//   move_copies(mypath)circle(r=.1,$fn=32);
 385// Example(2D): Requested point total is 14 but 15 points output due to extra end point
 386//   mypath = subdivide_path(pentagon(side=2), [3,4,3,4], method="segment", closed=false);
 387//   move_copies(mypath)circle(r=.1,$fn=32);
 388// Example(2D): Since 17 is not divisible by 5, a completely uniform distribution is not possible. 
 389//   mypath = subdivide_path(pentagon(side=2), 17);
 390//   move_copies(mypath)circle(r=.1,$fn=32);
 391// Example(2D): With `exact=false` a uniform distribution, but only 15 points
 392//   mypath = subdivide_path(pentagon(side=2), 17, exact=false);
 393//   move_copies(mypath)circle(r=.1,$fn=32);
 394// Example(2D): With `exact=false` you can also get extra points, here 20 instead of requested 18
 395//   mypath = subdivide_path(pentagon(side=2), 18, exact=false);
 396//   move_copies(mypath)circle(r=.1,$fn=32);
 397// Example(2D): Using refine in this example multiplies the point count by 3 by adding 2 points to each edge
 398//   mypath = subdivide_path(pentagon(side=2), refine=3);
 399//   move_copies(mypath)circle(r=.1,$fn=32);
 400// Example(2D): But note that refine doesn't distribute evenly by segment unless you change the method.  with the default method set to `"length"`, the points are distributed with more on the long segments in this example using refine.  
 401//   mypath = subdivide_path(square([8,2],center=true), refine=3);
 402//   move_copies(mypath)circle(r=.2,$fn=32);
 403// Example(2D): In this example with maxlen, every side gets a different number of new points
 404//   path = [[0,0],[0,4],[10,6],[10,0]];
 405//   spath = subdivide_path(path, maxlen=2, closed=true);
 406//   move_copies(spath) circle(r=.25,$fn=12);
 407// Example(FlatSpin,VPD=15,VPT=[0,0,1.5]): Three-dimensional paths also work
 408//   mypath = subdivide_path([[0,0,0],[2,0,1],[2,3,2]], 12);
 409//   move_copies(mypath)sphere(r=.1,$fn=32);
 410function subdivide_path(path, n, refine, maxlen, closed=true, exact, method) =
 411    let(path = force_path(path))
 412    assert(is_path(path))
 413    assert(num_defined([n,refine,maxlen]),"Must give exactly one of n, refine, and maxlen")
 414    refine==1 || n==len(path) ? path :
 415    is_def(maxlen) ?
 416        assert(is_undef(method), "Cannot give method with maxlen")
 417        assert(is_undef(exact), "Cannot give exact with maxlen")
 418        [
 419         for (p=pair(path,closed))
 420           let(steps = ceil(norm(p[1]-p[0])/maxlen))
 421           each lerpn(p[0], p[1], steps, false),
 422         if (!closed) last(path)
 423        ]               
 424    :
 425    let(
 426        exact = default(exact, true),
 427        method = default(method, "length")
 428    )
 429    assert(method=="length" || method=="segment")
 430    let(
 431        n = !is_undef(n)? n :
 432            !is_undef(refine)? len(path) * refine :
 433            undef
 434    )
 435    assert((is_num(n) && n>0) || is_vector(n),"Parameter n to subdivide_path must be postive number or vector")
 436    let(
 437        count = len(path) - (closed?0:1), 
 438        add_guess = method=="segment"?
 439                       (
 440                          is_list(n)
 441                          ? assert(len(n)==count,"Vector parameter n to subdivide_path has the wrong length")
 442                            add_scalar(n,-1)
 443                          : repeat((n-len(path)) / count, count)
 444                       )
 445                  : // method=="length"
 446                    assert(is_num(n),"Parameter n to subdivide path must be a number when method=\"length\"")
 447                    let(
 448                        path_lens = path_segment_lengths(path,closed),
 449                        add_density = (n - len(path)) / sum(path_lens)
 450                    )
 451                    path_lens * add_density,
 452        add = exact? _sum_preserving_round(add_guess)
 453                   : [for (val=add_guess) round(val)]
 454    )
 455    [
 456        for (i=[0:1:count-1]) 
 457           each lerpn(path[i],select(path,i+1), 1+add[i],endpoint=false),
 458        if (!closed) last(path)
 459    ];
 460
 461
 462
 463
 464// Function: resample_path()
 465// Synopsis: Returns an equidistant set of points along a path.
 466// SynTags: Path
 467// Topics: Paths
 468// See Also: subdivide_path()
 469// Usage:
 470//   newpath = resample_path(path, n|spacing=, [closed=]);
 471// Description:
 472//   Compute a uniform resampling of the input path.  If you specify `n` then the output path will have n
 473//   points spaced uniformly (by linear interpolation along the input path segments).  The only points of the
 474//   input path that are guaranteed to appear in the output path are the starting and ending points, and any
 475//   points that have an angular deflection of at least the number of degrees given in `keep_corners`.
 476//   If you specify `spacing` then the length you give will be rounded to the nearest spacing that gives
 477//   a uniform sampling of the path and the resulting uniformly sampled path is returned.
 478//   Note that because this function operates on a discrete input path the quality of the output depends on
 479//   the sampling of the input.  If you want very accurate output, use a lot of points for the input.
 480// Arguments:
 481//   path = path in any dimension or a 1-region
 482//   n = Number of points in output
 483//   ---
 484//   spacing = Approximate spacing desired
 485//   keep_corners = If given a scalar, path vertices with deflection angle greater than this are preserved in the output.
 486//   closed = set to true if path is closed.  Default: true
 487// Example(2D):  Subsampling lots of points from a smooth curve
 488//   path = xscale(2,circle($fn=250, r=10));
 489//   sampled = resample_path(path, 16);
 490//   stroke(path);
 491//   color("red")move_copies(sampled) circle($fn=16);
 492// Example(2D): Specified spacing is rounded to make a uniform sampling
 493//   path = xscale(2,circle($fn=250, r=10));
 494//   sampled = resample_path(path, spacing=17);
 495//   stroke(path);
 496//   color("red")move_copies(sampled) circle($fn=16);
 497// Example(2D): Notice that the corners are excluded.
 498//   path = square(20);
 499//   sampled = resample_path(path, spacing=6);
 500//   stroke(path,closed=true);
 501//   color("red")move_copies(sampled) circle($fn=16);
 502// Example(2D): Forcing preservation of corners.
 503//   path = square(20);
 504//   sampled = resample_path(path, spacing=6, keep_corners=90);
 505//   stroke(path,closed=true);
 506//   color("red")move_copies(sampled) circle($fn=16);
 507// Example(2D): Closed set to false
 508//   path = square(20);
 509//   sampled = resample_path(path, spacing=6,closed=false);
 510//   stroke(path);
 511//   color("red")move_copies(sampled) circle($fn=16);
 512
 513function resample_path(path, n, spacing, keep_corners, closed=true) =
 514    let(path = force_path(path))
 515    assert(is_path(path))
 516    assert(num_defined([n,spacing])==1,"Must define exactly one of n and spacing")
 517    assert(n==undef || (is_integer(n) && n>0))
 518    assert(spacing==undef || (is_finite(spacing) && spacing>0))
 519    assert(is_bool(closed))
 520    let(
 521        corners = is_undef(keep_corners)
 522          ? [0, len(path)-(closed?0:1)]
 523          : [
 524                0,
 525                for (i = [1:1:len(path)-(closed?1:2)])
 526                    let( ang = abs(modang(vector_angle(select(path,i-1,i+1))-180)) )
 527                    if (ang >= keep_corners) i,
 528                len(path)-(closed?0:1),
 529            ],
 530        pcnt = len(path),
 531        plen = path_length(path, closed=closed),
 532        subpaths = [ for (p = pair(corners)) [for(i = [p.x:1:p.y]) path[i%pcnt]] ],
 533        n = is_undef(n)? undef : closed? n+1 : n
 534    )
 535    assert(n==undef || n >= len(corners), "There are nore than `n=` corners whose angle is greater than `keep_corners=`.")
 536    let(
 537        lens = [for (subpath = subpaths) path_length(subpath)],
 538        part_ns = is_undef(n)
 539          ? [for (i=idx(subpaths)) max(1,round(lens[i]/spacing)-1)]
 540          : let(
 541                ccnt = len(corners),
 542                parts = [for (l=lens) (n-ccnt) * l/plen]
 543            )
 544            _sum_preserving_round(parts),
 545        out = [
 546            for (i = idx(subpaths))
 547                let(
 548                    subpath = subpaths[i],
 549                    splen = lens[i],
 550                    pn = part_ns[i] + 1,
 551                    distlist = lerpn(0, splen, pn, false),
 552                    cuts = path_cut_points(subpath, distlist, closed=false)
 553                )
 554                each column(cuts,0),
 555            if (!closed) last(path)
 556        ]
 557    ) out;
 558
 559
 560// Section: Path Geometry
 561
 562// Function: is_path_simple()
 563// Synopsis: Returns true if a path has no self intersections.
 564// Topics: Paths
 565// See Also: is_path()
 566// Usage:
 567//   bool = is_path_simple(path, [closed], [eps]);
 568// Description:
 569//   Returns true if the given 2D path is simple, meaning that it has no self-intersections.
 570//   Repeated points are not considered self-intersections: a path with such points can
 571//   still be simple.  
 572//   If closed is set to true then treat the path as a polygon.
 573// Arguments:
 574//   path = 2D path or 1-region
 575//   closed = set to true to treat path as a polygon.  Default: false
 576//   eps = Epsilon error value used for determine if points coincide.  Default: `EPSILON` (1e-9)
 577function is_path_simple(path, closed, eps=EPSILON) =
 578    is_1region(path) ? is_path_simple(path[0], default(closed,true), eps) :
 579    let(closed=default(closed,false))
 580    assert(is_path(path, 2),"Must give a 2D path")
 581    assert(is_bool(closed))
 582    let(
 583        path = deduplicate(path,closed=closed,eps=eps)
 584    )
 585    // check for path reversals
 586    [for(i=[0:1:len(path)-(closed?2:3)])
 587         let(v1=path[i+1]-path[i],
 588             v2=select(path,i+2)-path[i+1],
 589             normv1 = norm(v1),
 590             normv2 = norm(v2)
 591             )
 592         if (approx(v1*v2/normv1/normv2,-1)) 1
 593    ]  == [] 
 594    &&
 595    _path_self_intersections(path,closed=closed,eps=eps) == [];
 596
 597
 598// Function: path_closest_point()
 599// Synopsis: Returns the closest place on a path to a given point.
 600// Topics: Paths
 601// See Also: point_line_distance(), line_closest_point()
 602// Usage:
 603//   index_pt = path_closest_point(path, pt);
 604// Description:
 605//   Finds the closest path segment, and point on that segment to the given point.
 606//   Returns `[SEGNUM, POINT]`
 607// Arguments:
 608//   path = Path of any dimension or a 1-region.
 609//   pt = The point to find the closest point to.
 610//   closed = If true, the path is considered closed.
 611// Example(2D):
 612//   path = circle(d=100,$fn=6);
 613//   pt = [20,10];
 614//   closest = path_closest_point(path, pt);
 615//   stroke(path, closed=true);
 616//   color("blue") translate(pt) circle(d=3, $fn=12);
 617//   color("red") translate(closest[1]) circle(d=3, $fn=12);
 618function path_closest_point(path, pt, closed=true) =
 619    let(path = force_path(path))
 620    assert(is_path(path), "Input must be a path")
 621    assert(is_vector(pt, len(path[0])), "Input pt must be a compatible vector")
 622    assert(is_bool(closed))
 623    let(
 624        pts = [for (seg=pair(path,closed)) line_closest_point(seg,pt,SEGMENT)],
 625        dists = [for (p=pts) norm(p-pt)],
 626        min_seg = min_index(dists)
 627    ) [min_seg, pts[min_seg]];
 628
 629
 630// Function: path_tangents()
 631// Synopsis: Returns tangent vectors for each point along a path.
 632// Topics: Paths
 633// See Also: path_normals()
 634// Usage:
 635//   tangs = path_tangents(path, [closed], [uniform]);
 636// Description:
 637//   Compute the tangent vector to the input path.  The derivative approximation is described in deriv().
 638//   The returns vectors will be normalized to length 1.  If any derivatives are zero then
 639//   the function fails with an error.  If you set `uniform` to false then the sampling is
 640//   assumed to be non-uniform and the derivative is computed with adjustments to produce corrected
 641//   values.
 642// Arguments:
 643//   path = path of any dimension or a 1-region
 644//   closed = set to true of the path is closed.  Default: false
 645//   uniform = set to false to correct for non-uniform sampling.  Default: true
 646// Example(2D): A shape with non-uniform sampling gives distorted derivatives that may be undesirable.  Note that derivatives tilt towards the long edges of the rectangle.  
 647//   rect = square([10,3]);
 648//   tangents = path_tangents(rect,closed=true);
 649//   stroke(rect,closed=true, width=0.25);
 650//   color("purple")
 651//       for(i=[0:len(tangents)-1])
 652//           stroke([rect[i]-tangents[i], rect[i]+tangents[i]],width=.25, endcap2="arrow2");
 653// Example(2D): Setting uniform to false corrects the distorted derivatives for this example:
 654//   rect = square([10,3]);
 655//   tangents = path_tangents(rect,closed=true,uniform=false);
 656//   stroke(rect,closed=true, width=0.25);
 657//   color("purple")
 658//       for(i=[0:len(tangents)-1])
 659//           stroke([rect[i]-tangents[i], rect[i]+tangents[i]],width=.25, endcap2="arrow2");
 660function path_tangents(path, closed, uniform=true) =
 661    is_1region(path) ? path_tangents(path[0], default(closed,true), uniform) :
 662    let(closed=default(closed,false))
 663    assert(is_bool(closed))
 664    assert(is_path(path))
 665    !uniform ? [for(t=deriv(path,closed=closed, h=path_segment_lengths(path,closed))) unit(t)]
 666             : [for(t=deriv(path,closed=closed)) unit(t)];
 667
 668
 669// Function: path_normals()
 670// Synopsis: Returns normal vectors for each point along a path.
 671// Topics: Paths
 672// See Also: path_tangents()
 673// Usage:
 674//   norms = path_normals(path, [tangents], [closed]);
 675// Description:
 676//   Compute the normal vector to the input path.  This vector is perpendicular to the
 677//   path tangent and lies in the plane of the curve.  For 3d paths we define the plane of the curve
 678//   at path point i to be the plane defined by point i and its two neighbors.  At the endpoints of open paths
 679//   we use the three end points.  For 3d paths the computed normal is the one lying in this plane that points
 680//   towards the center of curvature at that path point.  For 2d paths, which lie in the xy plane, the normal
 681//   is the path pointing to the right of the direction the path is traveling.  If points are collinear then
 682//   a 3d path has no center of curvature, and hence the 
 683//   normal is not uniquely defined.  In this case the function issues an error.
 684//   For 2d paths the plane is always defined so the normal fails to exist only
 685//   when the derivative is zero (in the case of repeated points).
 686// Arguments:
 687//   path = 2D or 3D path or a 1-region
 688//   tangents = path tangents optionally supplied
 689//   closed = if true path is treated as a polygon.  Default: false
 690function path_normals(path, tangents, closed) =
 691    is_1region(path) ? path_normals(path[0], tangents, default(closed,true)) :
 692    let(closed=default(closed,false))
 693    assert(is_path(path,[2,3]))
 694    assert(is_bool(closed))
 695    let(
 696         tangents = default(tangents, path_tangents(path,closed)),
 697         dim=len(path[0])
 698    )
 699    assert(is_path(tangents) && len(tangents[0])==dim,"Dimensions of path and tangents must match")
 700    [
 701     for(i=idx(path))
 702         let(
 703             pts = i==0 ? (closed? select(path,-1,1) : select(path,0,2))
 704                 : i==len(path)-1 ? (closed? select(path,i-1,i+1) : select(path,i-2,i))
 705                 : select(path,i-1,i+1)
 706        )
 707        dim == 2 ? [tangents[i].y,-tangents[i].x]
 708                 : let( v=cross(cross(pts[1]-pts[0], pts[2]-pts[0]),tangents[i]))
 709                   assert(norm(v)>EPSILON, "3D path contains collinear points")
 710                   unit(v)
 711    ];
 712
 713
 714// Function: path_curvature()
 715// Synopsis: Returns the estimated numerical curvature of the path.
 716// Topics: Paths
 717// See Also: path_tangents(), path_normals(), path_torsion()
 718// Usage:
 719//   curvs = path_curvature(path, [closed]);
 720// Description:
 721//   Numerically estimate the curvature of the path (in any dimension).
 722// Arguments:
 723//   path = path in any dimension or a 1-region
 724//   closed = if true then treat the path as a polygon.  Default: false
 725function path_curvature(path, closed) =
 726    is_1region(path) ? path_curvature(path[0], default(closed,true)) :
 727    let(closed=default(closed,false))
 728    assert(is_bool(closed))
 729    assert(is_path(path))
 730    let( 
 731        d1 = deriv(path, closed=closed),
 732        d2 = deriv2(path, closed=closed)
 733    ) [
 734        for(i=idx(path))
 735        sqrt(
 736            sqr(norm(d1[i])*norm(d2[i])) -
 737            sqr(d1[i]*d2[i])
 738        ) / pow(norm(d1[i]),3)
 739    ];
 740
 741
 742// Function: path_torsion()
 743// Synopsis: Returns the estimated numerical torsion of the path.
 744// Topics: Paths
 745// See Also: path_tangents(), path_normals(), path_curvature()
 746// Usage:
 747//   torsions = path_torsion(path, [closed]);
 748// Description:
 749//   Numerically estimate the torsion of a 3d path.
 750// Arguments:
 751//   path = 3D path
 752//   closed = if true then treat path as a polygon.  Default: false
 753function path_torsion(path, closed=false) =
 754    assert(is_path(path,3), "Input path must be a 3d path")
 755    assert(is_bool(closed))
 756    let(
 757        d1 = deriv(path,closed=closed),
 758        d2 = deriv2(path,closed=closed),
 759        d3 = deriv3(path,closed=closed)
 760    ) [
 761        for (i=idx(path)) let(
 762            crossterm = cross(d1[i],d2[i])
 763        ) crossterm * d3[i] / sqr(norm(crossterm))
 764    ];
 765
 766
 767// Section: Breaking paths up into subpaths
 768
 769
 770
 771// Function: path_cut()
 772// Synopsis: Cuts a path into subpaths at various points.
 773// SynTags: PathList
 774// Topics: Paths, Path Subdivision
 775// See Also: split_path_at_self_crossings(), path_cut_points()
 776// Usage:
 777//   path_list = path_cut(path, cutdist, [closed]);
 778// Description:
 779//   Given a list of distances in `cutdist`, cut the path into
 780//   subpaths at those lengths, returning a list of paths.
 781//   If the input path is closed then the final path will include the
 782//   original starting point.  The list of cut distances must be
 783//   in ascending order and should not include the endpoints: 0 
 784//   or len(path).  If you repeat a distance you will get an
 785//   empty list in that position in the output.  If you give an
 786//   empty cutdist array you will get the input path as output
 787//   (without the final vertex doubled in the case of a closed path).
 788// Arguments:
 789//   path = path of any dimension or a 1-region
 790//   cutdist = Distance or list of distances where path is cut
 791//   closed = If true, treat the path as a closed polygon.  Default: false
 792// Example(2D,NoAxes):
 793//   path = circle(d=100);
 794//   segs = path_cut(path, [50, 200], closed=true);
 795//   rainbow(segs) stroke($item, endcaps="butt", width=3);
 796function path_cut(path,cutdist,closed) =
 797  is_num(cutdist) ? path_cut(path,[cutdist],closed) :
 798  is_1region(path) ? path_cut(path[0], cutdist, default(closed,true)):
 799  let(closed=default(closed,false))
 800  assert(is_bool(closed))
 801  assert(is_vector(cutdist))
 802  assert(last(cutdist)<path_length(path,closed=closed)-EPSILON,"Cut distances must be smaller than the path length")
 803  assert(cutdist[0]>EPSILON, "Cut distances must be strictly positive")
 804  let(
 805      cutlist = path_cut_points(path,cutdist,closed=closed)
 806  )
 807  _path_cut_getpaths(path, cutlist, closed);
 808
 809
 810function _path_cut_getpaths(path, cutlist, closed) =
 811  let(
 812      cuts = len(cutlist)
 813  )
 814  [
 815      [ each list_head(path,cutlist[0][1]-1),
 816        if (!approx(cutlist[0][0], path[cutlist[0][1]-1])) cutlist[0][0]
 817      ],
 818      for(i=[0:1:cuts-2])
 819          cutlist[i][0]==cutlist[i+1][0] && cutlist[i][1]==cutlist[i+1][1] ? []
 820          :
 821          [ if (!approx(cutlist[i][0], select(path,cutlist[i][1]))) cutlist[i][0],
 822            each slice(path, cutlist[i][1], cutlist[i+1][1]-1),
 823            if (!approx(cutlist[i+1][0], select(path,cutlist[i+1][1]-1))) cutlist[i+1][0],
 824          ],
 825      [
 826        if (!approx(cutlist[cuts-1][0], select(path,cutlist[cuts-1][1]))) cutlist[cuts-1][0],
 827        each select(path,cutlist[cuts-1][1],closed ? 0 : -1)
 828      ]
 829  ];
 830
 831
 832
 833// Function: path_cut_points()
 834// Synopsis: Returns a list of cut points at a list of distances from the first point in a path.
 835// Topics: Paths, Path Subdivision
 836// See Also: path_cut(), split_path_at_self_crossings()
 837// Usage:
 838//   cuts = path_cut_points(path, cutdist, [closed=], [direction=]);
 839//
 840// Description:
 841//   Cuts a path at a list of distances from the first point in the path.  Returns a list of the cut
 842//   points and indices of the next point in the path after that point.  So for example, a return
 843//   value entry of [[2,3], 5] means that the cut point was [2,3] and the next point on the path after
 844//   this point is path[5].  If the path is too short then path_cut_points returns undef.  If you set
 845//   `direction` to true then `path_cut_points` will also return the tangent vector to the path and a normal
 846//   vector to the path.  It tries to find a normal vector that is coplanar to the path near the cut
 847//   point.  If this fails it will return a normal vector parallel to the xy plane.  The output with
 848//   direction vectors will be `[point, next_index, tangent, normal]`.
 849//   .
 850//   If you give the very last point of the path as a cut point then the returned index will be
 851//   one larger than the last index (so it will not be a valid index).  If you use the closed
 852//   option then the returned index will be equal to the path length for cuts along the closing
 853//   path segment, and if you give a point equal to the path length you will get an
 854//   index of len(path)+1 for the index.  
 855//
 856// Arguments:
 857//   path = path to cut
 858//   cutdist = distances where the path should be cut (a list) or a scalar single distance
 859//   ---
 860//   closed = set to true if the curve is closed.  Default: false
 861//   direction = set to true to return direction vectors.  Default: false
 862//
 863// Example(NORENDER):
 864//   square=[[0,0],[1,0],[1,1],[0,1]];
 865//   path_cut_points(square, [.5,1.5,2.5]);   // Returns [[[0.5, 0], 1], [[1, 0.5], 2], [[0.5, 1], 3]]
 866//   path_cut_points(square, [0,1,2,3]);      // Returns [[[0, 0], 1], [[1, 0], 2], [[1, 1], 3], [[0, 1], 4]]
 867//   path_cut_points(square, [0,0.8,1.6,2.4,3.2], closed=true);  // Returns [[[0, 0], 1], [[0.8, 0], 1], [[1, 0.6], 2], [[0.6, 1], 3], [[0, 0.8], 4]]
 868//   path_cut_points(square, [0,0.8,1.6,2.4,3.2]);               // Returns [[[0, 0], 1], [[0.8, 0], 1], [[1, 0.6], 2], [[0.6, 1], 3], undef]
 869function path_cut_points(path, cutdist, closed=false, direction=false) =
 870    let(long_enough = len(path) >= (closed ? 3 : 2))
 871    assert(long_enough,len(path)<2 ? "Two points needed to define a path" : "Closed path must include three points")
 872    is_num(cutdist) ? path_cut_points(path, [cutdist],closed, direction)[0] :
 873    assert(is_vector(cutdist))
 874    assert(is_increasing(cutdist), "Cut distances must be an increasing list")
 875    let(cuts = path_cut_points_recurse(path,cutdist,closed))
 876    !direction
 877       ? cuts
 878       : let(
 879             dir = _path_cuts_dir(path, cuts, closed),
 880             normals = _path_cuts_normals(path, cuts, dir, closed)
 881         )
 882         hstack(cuts, list_to_matrix(dir,1), list_to_matrix(normals,1));
 883
 884// Main recursive path cut function
 885function path_cut_points_recurse(path, dists, closed=false, pind=0, dtotal=0, dind=0, result=[]) =
 886    dind == len(dists) ? result :
 887    let(
 888        lastpt = len(result)==0? [] : last(result)[0],       // location of last cut point
 889        dpartial = len(result)==0? 0 : norm(lastpt-select(path,pind)),  // remaining length in segment
 890        nextpoint = dists[dind] < dpartial+dtotal  // Do we have enough length left on the current segment?
 891           ? [lerp(lastpt,select(path,pind),(dists[dind]-dtotal)/dpartial),pind] 
 892           : _path_cut_single(path, dists[dind]-dtotal-dpartial, closed, pind)
 893    ) 
 894    path_cut_points_recurse(path, dists, closed, nextpoint[1], dists[dind],dind+1, concat(result, [nextpoint]));
 895
 896
 897// Search for a single cut point in the path
 898function _path_cut_single(path, dist, closed=false, ind=0, eps=1e-7) =
 899    // If we get to the very end of the path (ind is last point or wraparound for closed case) then
 900    // check if we are within epsilon of the final path point.  If not we're out of path, so we fail
 901    ind==len(path)-(closed?0:1) ?
 902       assert(dist<eps,"Path is too short for specified cut distance")
 903       [select(path,ind),ind+1]
 904    :let(d = norm(path[ind]-select(path,ind+1))) d > dist ?
 905        [lerp(path[ind],select(path,ind+1),dist/d), ind+1] :
 906        _path_cut_single(path, dist-d,closed, ind+1, eps);
 907
 908// Find normal directions to the path, coplanar to local part of the path
 909// Or return a vector parallel to the x-y plane if the above fails
 910function _path_cuts_normals(path, cuts, dirs, closed=false) =
 911    [for(i=[0:len(cuts)-1])
 912        len(path[0])==2? [-dirs[i].y, dirs[i].x]
 913          : 
 914            let(
 915                plane = len(path)<3 ? undef :
 916                let(start = max(min(cuts[i][1],len(path)-1),2)) _path_plane(path, start, start-2)
 917            )
 918            plane==undef?
 919                ( dirs[i].x==0 && dirs[i].y==0 ? [1,0,0]  // If it's z direction return x vector
 920                                               : unit([-dirs[i].y, dirs[i].x,0])) // otherwise perpendicular to projection
 921                : unit(cross(dirs[i],cross(plane[0],plane[1])))
 922    ];
 923
 924// Scan from the specified point (ind) to find a noncoplanar triple to use
 925// to define the plane of the path.
 926function _path_plane(path, ind, i,closed) =
 927    i<(closed?-1:0) ? undef :
 928    !is_collinear(path[ind],path[ind-1], select(path,i))?
 929        [select(path,i)-path[ind-1],path[ind]-path[ind-1]] :
 930        _path_plane(path, ind, i-1);
 931
 932// Find the direction of the path at the cut points
 933function _path_cuts_dir(path, cuts, closed=false, eps=1e-2) =
 934    [for(ind=[0:len(cuts)-1])
 935        let(
 936            zeros = path[0]*0,
 937            nextind = cuts[ind][1],
 938            nextpath = unit(select(path, nextind+1)-select(path, nextind),zeros),
 939            thispath = unit(select(path, nextind) - select(path,nextind-1),zeros),
 940            lastpath = unit(select(path,nextind-1) - select(path, nextind-2),zeros),
 941            nextdir =
 942                nextind==len(path) && !closed? lastpath :
 943                (nextind<=len(path)-2 || closed) && approx(cuts[ind][0], path[nextind],eps)
 944                   ? unit(nextpath+thispath)
 945              : (nextind>1 || closed) && approx(cuts[ind][0],select(path,nextind-1),eps)
 946                   ? unit(thispath+lastpath)
 947              :  thispath
 948        ) nextdir
 949    ];
 950
 951
 952// internal function
 953// converts pathcut output form to a [segment, u]
 954// form list that works withi path_select
 955function _cut_to_seg_u_form(pathcut, path, closed) =
 956  let(lastind = len(path) - (closed?0:1))
 957  [for(entry=pathcut)
 958    entry[1] > lastind ? [lastind,0] :
 959    let(
 960        a = path[entry[1]-1],
 961        b = path[entry[1]],
 962        c = entry[0],
 963        i = max_index(v_abs(b-a)),
 964        factor = (c[i]-a[i])/(b[i]-a[i])
 965    )
 966    [entry[1]-1,factor]
 967  ];
 968
 969
 970
 971// Function: split_path_at_self_crossings()
 972// Synopsis: Split a 2D path wherever it crosses itself.
 973// SynTags: PathList
 974// Topics: Paths, Path Subdivision
 975// See Also: path_cut(), path_cut_points()
 976// Usage:
 977//   paths = split_path_at_self_crossings(path, [closed], [eps]);
 978// Description:
 979//   Splits a 2D path into sub-paths wherever the original path crosses itself.
 980//   Splits may occur mid-segment, so new vertices will be created at the intersection points.
 981//   Returns a list of the resulting subpaths.  
 982// Arguments:
 983//   path = A 2D path or a 1-region.
 984//   closed = If true, treat path as a closed polygon.  Default: true
 985//   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 986// Example(2D,NoAxes):
 987//   path = [ [-100,100], [0,-50], [100,100], [100,-100], [0,50], [-100,-100] ];
 988//   paths = split_path_at_self_crossings(path);
 989//   rainbow(paths) stroke($item, closed=false, width=3);
 990function split_path_at_self_crossings(path, closed=true, eps=EPSILON) =
 991    let(path = force_path(path))
 992    assert(is_path(path,2), "Must give a 2D path")
 993    assert(is_bool(closed))
 994    let(
 995        path = list_unwrap(path, eps=eps),
 996        isects = deduplicate(
 997            eps=eps,
 998            concat(
 999                [[0, 0]],
1000                sort([
1001                    for (
1002                        a = _path_self_intersections(path, closed=closed, eps=eps),
1003                        ss = [ [a[1],a[2]], [a[3],a[4]] ]
1004                    ) if (ss[0] != undef) ss
1005                ]),
1006                [[len(path)-(closed?1:2), 1]]
1007            )
1008        )
1009    ) [
1010        for (p = pair(isects))
1011            let(
1012                s1 = p[0][0],
1013                u1 = p[0][1],
1014                s2 = p[1][0],
1015                u2 = p[1][1],
1016                section = _path_select(path, s1, u1, s2, u2, closed=closed),
1017                outpath = deduplicate(eps=eps, section)
1018            )
1019            if (len(outpath)>1) outpath
1020    ];
1021
1022
1023function _tag_self_crossing_subpaths(path, nonzero, closed=true, eps=EPSILON) =
1024    let(
1025        subpaths = split_path_at_self_crossings(
1026            path, closed=true, eps=eps
1027        )
1028    ) [
1029        for (subpath = subpaths) let(
1030            seg = select(subpath,0,1),
1031            mp = mean(seg),
1032            n = line_normal(seg) / 2048,
1033            p1 = mp + n,
1034            p2 = mp - n,
1035            p1in = point_in_polygon(p1, path, nonzero=nonzero) >= 0,
1036            p2in = point_in_polygon(p2, path, nonzero=nonzero) >= 0,
1037            tag = (p1in && p2in)? "I" : "O"
1038        ) [tag, subpath]
1039    ];
1040
1041
1042// Function: polygon_parts()
1043// Synopsis: Parses a self-intersecting polygon into a list of non-intersecting polygons.
1044// SynTags: PathList
1045// Topics: Paths, Polygons
1046// See Also: split_path_at_self_crossings(), path_cut(), path_cut_points()
1047// Usage:
1048//   splitpolys = polygon_parts(poly, [nonzero], [eps]);
1049// Description:
1050//   Given a possibly self-intersecting 2d polygon, constructs a representation of the original polygon as a list of
1051//   non-intersecting simple polygons.  If nonzero is set to true then it uses the nonzero method for defining polygon membership.
1052//   For simple cases, such as the pentagram, this will produce the outer perimeter of a self-intersecting polygon.  
1053// Arguments:
1054//   poly = a 2D polygon or 1-region
1055//   nonzero = If true use the nonzero method for checking if a point is in a polygon.  Otherwise use the even-odd method.  Default: false
1056//   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
1057// Example(2D,NoAxes):  This cross-crossing polygon breaks up into its 3 components (regardless of the value of nonzero).
1058//   poly = [
1059//       [-100,100], [0,-50], [100,100],
1060//       [100,-100], [0,50], [-100,-100]
1061//   ];
1062//   splitpolys = polygon_parts(poly);
1063//   rainbow(splitpolys) stroke($item, closed=true, width=3);
1064// Example(2D,NoAxes): With nonzero=false you get even-odd mode which matches OpenSCAD, so the pentagram breaks apart into its five points.
1065//   pentagram = turtle(["move",100,"left",144], repeat=4);
1066//   left(100)polygon(pentagram);
1067//   rainbow(polygon_parts(pentagram,nonzero=false))
1068//     stroke($item,closed=true,width=2.5);
1069// Example(2D,NoAxes): With nonzero=true you get only the outer perimeter.  You can use this to create the polygon using the nonzero method, which is not supported by OpenSCAD.
1070//   pentagram = turtle(["move",100,"left",144], repeat=4);
1071//   outside = polygon_parts(pentagram,nonzero=true);
1072//   left(100)region(outside);
1073//   rainbow(outside)
1074//     stroke($item,closed=true,width=2.5);
1075// Example(2D,NoAxes): 
1076//   N=12;
1077//   ang=360/N;
1078//   sr=10;
1079//   poly = turtle(["angle", 90+ang/2,
1080//                  "move", sr, "left",
1081//                  "move", 2*sr*sin(ang/2), "left",
1082//                  "repeat", 4,
1083//                     ["move", 2*sr, "left",
1084//                      "move", 2*sr*sin(ang/2), "left"],
1085//                  "move", sr]);
1086//   stroke(poly, width=.3);
1087//   right(20)rainbow(polygon_parts(poly)) polygon($item);
1088// Example(2D,NoAxes): overlapping poly segments disappear
1089//   poly = [[0,0], [10,0], [10,10], [0,10],[0,20], [20,10],[10,10], [0,10],[0,0]];
1090//   stroke(poly,width=0.3);
1091//   right(22)stroke(polygon_parts(poly)[0], width=0.3, closed=true);
1092// Example(2D,NoAxes): Poly segments disappear outside as well
1093//   poly = turtle(["repeat", 3, ["move", 17, "left", "move", 10, "left", "move", 7, "left", "move", 10, "left"]]);
1094//   back(2)stroke(poly,width=.5);
1095//   fwd(12)rainbow(polygon_parts(poly)) stroke($item, closed=true, width=0.5);
1096// Example(2D,NoAxes):  This shape has six components
1097//   poly = turtle(["repeat", 3, ["move", 15, "left", "move", 7, "left", "move", 10, "left", "move", 17, "left"]]);
1098//   polygon(poly);
1099//   right(22)rainbow(polygon_parts(poly)) polygon($item);
1100// Example(2D,NoAxes): When the loops of the shape overlap then nonzero gives a different result than the even-odd method.
1101//   poly = turtle(["repeat", 3, ["move", 15, "left", "move", 7, "left", "move", 10, "left", "move", 10, "left"]]);
1102//   polygon(poly);
1103//   right(27)rainbow(polygon_parts(poly)) polygon($item);
1104//   move([16,-14])rainbow(polygon_parts(poly,nonzero=true)) polygon($item);
1105function polygon_parts(poly, nonzero=false, eps=EPSILON) =
1106    let(poly = force_path(poly))
1107    assert(is_path(poly,2), "Must give 2D polygon")
1108    assert(is_bool(nonzero))    
1109    let(
1110        poly = list_unwrap(poly, eps=eps),
1111        tagged = _tag_self_crossing_subpaths(poly, nonzero=nonzero, closed=true, eps=eps),
1112        kept = [for (sub = tagged) if(sub[0] == "O") sub[1]],
1113        outregion = _assemble_path_fragments(kept, eps=eps)
1114    ) outregion;
1115
1116
1117function _extreme_angle_fragment(seg, fragments, rightmost=true, eps=EPSILON) =
1118    !fragments? [undef, []] :
1119    let(
1120        delta = seg[1] - seg[0],
1121        segang = atan2(delta.y,delta.x),
1122        frags = [
1123            for (i = idx(fragments)) let(
1124                fragment = fragments[i],
1125                fwdmatch = approx(seg[1], fragment[0], eps=eps),
1126                bakmatch =  approx(seg[1], last(fragment), eps=eps)
1127            ) [
1128                fwdmatch,
1129                bakmatch,
1130                bakmatch? reverse(fragment) : fragment
1131            ]
1132        ],
1133        angs = [
1134            for (frag = frags)
1135                (frag[0] || frag[1])? let(
1136                    delta2 = frag[2][1] - frag[2][0],
1137                    segang2 = atan2(delta2.y, delta2.x)
1138                ) modang(segang2 - segang) : (
1139                    rightmost? 999 : -999
1140                )
1141        ],
1142        fi = rightmost? min_index(angs) : max_index(angs)
1143    ) abs(angs[fi]) > 360? [undef, fragments] : let(
1144        remainder = [for (i=idx(fragments)) if (i!=fi) fragments[i]],
1145        frag = frags[fi],
1146        foundfrag = frag[2]
1147    ) [foundfrag, remainder];
1148
1149
1150/// Internal Function: _assemble_a_path_from_fragments()
1151/// Usage:
1152///   _assemble_a_path_from_fragments(subpaths);
1153/// Description:
1154///   Given a list of paths, assembles them together into one complete closed polygon path, and
1155///   remainder fragments.  Returns [PATH, FRAGMENTS] where FRAGMENTS is the list of remaining
1156///   unused path fragments.
1157/// Arguments:
1158///   fragments = List of paths to be assembled into complete polygons.
1159///   rightmost = If true, assemble paths using rightmost turns. Leftmost if false.
1160///   startfrag = The fragment to start with.  Default: 0
1161///   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
1162function _assemble_a_path_from_fragments(fragments, rightmost=true, startfrag=0, eps=EPSILON) =
1163    len(fragments)==0? [[],[]] :
1164    len(fragments)==1? [fragments[0],[]] :
1165    let(
1166        path = fragments[startfrag],
1167        newfrags = [for (i=idx(fragments)) if (i!=startfrag) fragments[i]]
1168    ) are_ends_equal(path, eps=eps)? (
1169        // starting fragment is already closed
1170        [path, newfrags]
1171    ) : let(
1172        // Find rightmost/leftmost continuation fragment
1173        seg = select(path,-2,-1),
1174        extrema = _extreme_angle_fragment(seg=seg, fragments=newfrags, rightmost=rightmost, eps=eps),
1175        foundfrag = extrema[0],
1176        remainder = extrema[1]
1177    ) is_undef(foundfrag)? (
1178        // No remaining fragments connect!  INCOMPLETE PATH!
1179        // Treat it as complete.
1180        [path, remainder]
1181    ) : are_ends_equal(foundfrag, eps=eps)? (
1182        // Found fragment is already closed
1183        [foundfrag, concat([path], remainder)]
1184    ) : let(
1185        fragend = last(foundfrag),
1186        hits = [for (i = idx(path,e=-2)) if(approx(path[i],fragend,eps=eps)) i]
1187    ) hits? (
1188        let(
1189            // Found fragment intersects with initial path
1190            hitidx = last(hits),
1191            newpath = list_head(path,hitidx),
1192            newfrags = concat(len(newpath)>1? [newpath] : [], remainder),
1193            outpath = concat(slice(path,hitidx,-2), foundfrag)
1194        )
1195        [outpath, newfrags]
1196    ) : let(
1197        // Path still incomplete.  Continue building it.
1198        newpath = concat(path, list_tail(foundfrag)),
1199        newfrags = concat([newpath], remainder)
1200    )
1201    _assemble_a_path_from_fragments(
1202        fragments=newfrags,
1203        rightmost=rightmost,
1204        eps=eps
1205    );
1206
1207
1208/// Internal Function: _assemble_path_fragments()
1209/// Usage:
1210///   _assemble_path_fragments(subpaths);
1211/// Description:
1212///   Given a list of paths, assembles them together into complete closed polygon paths if it can.
1213///   Polygons with area < eps will be discarded and not returned.  
1214/// Arguments:
1215///   fragments = List of paths to be assembled into complete polygons.
1216///   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
1217function _assemble_path_fragments(fragments, eps=EPSILON, _finished=[]) =
1218    len(fragments)==0? _finished :
1219    let(
1220        minxidx = min_index([
1221            for (frag=fragments) min(column(frag,0))
1222        ]),
1223        result_l = _assemble_a_path_from_fragments(
1224            fragments=fragments,
1225            startfrag=minxidx,
1226            rightmost=false,
1227            eps=eps
1228        ),
1229        result_r = _assemble_a_path_from_fragments(
1230            fragments=fragments,
1231            startfrag=minxidx,
1232            rightmost=true,
1233            eps=eps
1234        ),
1235        l_area = abs(polygon_area(result_l[0])),
1236        r_area = abs(polygon_area(result_r[0])),
1237        result = l_area < r_area? result_l : result_r,
1238        newpath = list_unwrap(result[0]),
1239        remainder = result[1],
1240        finished = min(l_area,r_area)<eps ? _finished : concat(_finished, [newpath])
1241    ) _assemble_path_fragments(
1242        fragments=remainder,
1243        eps=eps,
1244        _finished=finished
1245    );
1246
1247
1248
1249// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap